| Example Program Blast Reports Parsing the output of BLAST call. 1 | #include <iostream>
| 2 | #include <seqan/blast.h>
| 3 |
| 4 | using namespace std;
| 5 | using namespace seqan;
| 6 |
| 7 |
| 8 | template <typename TFile>
| 9 | void read_blast_report(TFile & strm)
| 10 | {
| 11 |
|
12 |
| 13 | typedef BlastHsp<BlastN, FullInfo> TBlastHsp;
| 14 |
|
15 |
| 16 | typedef BlastReport<TBlastHsp, StreamReport<TFile> > TBlastReport;
| 17 |
|
18 | typedef typename Hit<TBlastReport>::Type TBlastHit;
| 19 | typedef typename Iterator<TBlastReport,HitIterator>::Type THitIterator;
| 20 | typedef typename Iterator<TBlastHit,HspIterator>::Type THspIterator;
| 21 |
| 22 |
|
23 | TBlastReport blast;
| 24 |
| 25 | // counters
| 26 | unsigned hspcount = 0;
| 27 | unsigned hitcount = 0;
| 28 | unsigned highsignif = 0;
| 29 |
| 30 | while(!atEnd(strm,blast))
| 31 | {
| 32 | /// get the current Blast report (there can be multiple reports in one file)
| 33 | read(strm,blast,Blast());
| 34 | std::cout << "Query: "<<getQueryName(blast) <<"\n";
| 35 | std::cout << "Database: "<<getDatabaseName(blast) <<"\n\n";
| 36 |
| 37 | /// iterate over hits
| 38 | THitIterator hit_it(blast);
| 39 | for(; !atEnd(strm,hit_it); goNext(strm,hit_it))
| 40 | {
| 41 | ++hitcount;
| 42 | TBlastHit hit = getValue(strm,hit_it);
| 43 | std::cout << " Hit: " <<name(hit) <<"\n\n";
| 44 |
| 45 | /// iterate over alignments (HSPs)
| 46 | THspIterator hsp_it(hit);
| 47 | for(; !atEnd(strm,hsp_it); goNext(strm,hsp_it))
| 48 | {
| 49 | ++hspcount;
| 50 | TBlastHsp hsp = getValue(strm,hsp_it);
| 51 |
| 52 | /// do something with the alignment, e.g.
| 53 | /// output score and length of alignment
| 54 | std::cout << " Score = " << getBitScore(hsp) << "\n";
| 55 | std::cout << " Length = " << length(hsp) << "\n\n";
| 56 | /// and count alignments with highly significant e-values
| 57 | if(getEValue(hsp)<0.01)
| 58 | ++highsignif;
| 59 |
| 60 | }
| 61 | }
| 62 | }
| 63 | std::cout <<"Total number of Hits: "<< hitcount<<std::endl;
| 64 | std::cout <<"Total number of HSPs: "<< hspcount<<std::endl;
| 65 | std::cout <<"Number of highly significant HSPs: "<< highsignif<<std::endl;
| 66 |
| 67 | }
| 68 |
| 69 |
| 70 | int main()
| 71 | {
| 72 |
| 73 | std::fstream strm;
| 74 | strm.open("ecoln.out",ios_base::in | ios_base::binary);
| 75 | read_blast_report(strm);
| 76 |
| 77 | return 0;
| 78 | }
|
|